Models of vegetation composition based on climate predictors
Based on code from Martin Holdrege’s Sagebrush-Fire paper
Alice Stears
2025-04-11
- Dependencies
- read in data
- Prep data
- Add a constant to the response variable (+1) so that models run…
- Identify the ecoregion and response variable type to use in this model run
- Subset the data to only include data for the ecoregion of interest
- Currently, subsampling data from the “Texas Coastal Plain”, since it’s quite different from other regions and is really messing with model fit
- Visualize the response variable
- Visualize the predictor variables
- Plot predictor vars against each other
- Predictor variables compared to binned response variables
- Standardize the predictor variables for the model-fitting process
- Model Fitting
- Visualizations of Model Predictions and Residuals – using best lambda model
- Cross-validation
- Using best lambda model
- Use terms from global model to re-fit and predict on different held out regions
- Figures show residuals for each of the models fit to held-out ecoregions
- Below are the RMSE and bias values for predictions made for each holdout level II ecoregion, compared to predictions from the global model for that same ecoregion
- Using 1se lambda model
- Use terms from global model to re-fit and predict on different held out regions
- Figures show residuals for each of the models fit to held-out ecoregions
- Below are the RMSE and bias values for predictions made for each holdout level II ecoregion, compared to predictions from the global model for that same ecoregion
- Using best lambda model
- Save output
- session info
The data consists of vegetation % cover by functional group from across CONUS (from AIM, FIA, LANDFIRE, and RAP), as well as climate variables from DayMet, which have been aggregated into mean interannual conditions accross multiple temporal windows.
Dependencies
User defined parameters
print(params)## $run
## [1] TRUE
##
## $test_run
## [1] FALSE
##
## $save_figs
## [1] FALSE
##
## $ecoregion
## [1] "shrubGrass"
##
## $response
## [1] "BareGroundCover"# set to true if want to run for a limited number of rows (i.e. for code testing)
test_run <- params$test_run
save_figs <- params$save_figs
response <- params$response
fit_sample <- TRUE # fit model to a sample of the data
n_train <- 5e4 # sample size of the training data
n_test <- 1e6 # sample size of the testing data (if this is too big the decile dotplot code throws memory errors)
run <- params$run# set option so resampled dataset created here reproduces earlier runs of this code with dplyr 1.0.10
source("../../../Functions/glmTransformsIterates.R")
source("../../../Functions/transformPreds.R")
source("../../../Functions/StepBeta_mine.R")
#source("src/fig_params.R")
#source("src/modeling_functions.R")
library(ggspatial)
library(terra)
library(tidyterra)
library(sf)
library(caret)
library(tidyverse)
library(GGally) # for ggpairs()
library(pdp) # for partial dependence plots
library(gridExtra)
library(knitr)
library(patchwork) # for figure insets etc.
library(ggtext)
library(StepBeta)
theme_set(theme_classic())
library(here)
library(rsample)
library(kableExtra)
library(glmnet)
library(USA.state.boundaries)read in data
Data compiled in the prepDataForModels.R script
here::i_am("Analysis/VegComposition/ModelFitting/02_ModelFitting.Rmd")
modDat <- readRDS( here("Data_processed", "CoverData", "DataForModels_spatiallyAveraged_withSoils_noSf.rds"))
## there are some values of the annual wet degree days 5th percentile that have -Inf?? change to lowest value for now?
modDat[is.infinite(modDat$annWetDegDays_5percentile_3yrAnom), "annWetDegDays_5percentile_3yrAnom"] <- -47.8
## same, but for annual water deficit 95th percentile
modDat[is.infinite(modDat$annWaterDeficit_95percentile_3yrAnom), "annWaterDeficit_95percentile_3yrAnom"] <- -600
# # Convert total cover variables into proportions (for later use in beta regression models) ... proportions are already scaled from zero to 1
# modDat <- modDat %>%
# mutate(TotalTreeCover = TotalTreeCover/100,
# CAMCover = CAMCover/100,
# TotalHerbaceousCover = TotalHerbaceousCover/100,
# BareGroundCover = BareGroundCover/100,
# ShrubCover = ShrubCover/100
# )
# For all response variables, make sure there are no 0s add or subtract .0001 from each, since the Gamma model framework can't handle that
modDat[modDat$TotalTreeCover == 0 & !is.na(modDat$TotalTreeCover), "TotalTreeCover"] <- 0.0001
modDat[modDat$CAMCover == 0 & !is.na(modDat$CAMCover), "CAMCover"] <- 0.0001
modDat[modDat$TotalHerbaceousCover == 0 & !is.na(modDat$TotalHerbaceousCover), "TotalHerbaceousCover"] <- 0.0001
modDat[modDat$BareGroundCover == 0 & !is.na(modDat$BareGroundCover), "BareGroundCover"] <- 0.0001
modDat[modDat$ShrubCover == 0 & !is.na(modDat$ShrubCover), "ShrubCover"] <- 0.0001
modDat[modDat$BroadleavedTreeCover_prop == 0 & !is.na(modDat$BroadleavedTreeCover_prop), "BroadleavedTreeCover_prop"] <- 0.0001
modDat[modDat$NeedleLeavedTreeCover_prop == 0 & !is.na(modDat$NeedleLeavedTreeCover_prop), "NeedleLeavedTreeCover_prop"] <- 0.0001
modDat[modDat$C4Cover_prop == 0 & !is.na(modDat$C4Cover_prop), "C4Cover_prop"] <- 0.0001
modDat[modDat$C3Cover_prop == 0 & !is.na(modDat$C3Cover_prop), "C3Cover_prop"] <- 0.0001
modDat[modDat$ForbCover_prop == 0 & !is.na(modDat$ForbCover_prop), "ForbCover_prop"] <- 0.0001
#
# modDat[modDat$TotalTreeCover ==1& !is.na(modDat$TotalTreeCover), "TotalTreeCover"] <- 0.999
# modDat[modDat$CAMCover ==1& !is.na(modDat$CAMCover), "CAMCover"] <- 0.999
# modDat[modDat$TotalHerbaceousCover ==1 & !is.na(modDat$TotalHerbaceousCover), "TotalHerbaceousCover"] <- 0.999
# modDat[modDat$BareGroundCover ==1& !is.na(modDat$BareGroundCover), "BareGroundCover"] <- 0.999
# modDat[modDat$ShrubCover ==1& !is.na(modDat$ShrubCover), "ShrubCover"] <- 0.999
# modDat[modDat$BroadleavedTreeCover_prop ==1& !is.na(modDat$BroadleavedTreeCover_prop), "BroadleavedTreeCover_prop"] <- 0.999
# modDat[modDat$NeedleLeavedTreeCover_prop ==1& !is.na(modDat$NeedleLeavedTreeCover_prop), "NeedleLeavedTreeCover_prop"] <- 0.999
# modDat[modDat$C4Cover_prop ==1& !is.na(modDat$C4Cover_prop), "C4Cover_prop"] <- 0.999
# modDat[modDat$C3Cover_prop ==1& !is.na(modDat$C3Cover_prop), "C3Cover_prop"] <- 0.999
# modDat[modDat$ForbCover_prop ==1& !is.na(modDat$ForbCover_prop), "ForbCover_prop"] <- 0.999Prep data
set.seed(1234)
modDat_1 <- modDat %>%
select(-c(prcp_annTotal:annVPD_min)) %>%
# mutate(Lon = st_coordinates(.)[,1],
# Lat = st_coordinates(.)[,2]) %>%
# st_drop_geometry() %>%
# filter(!is.na(newRegion))
rename("tmin" = tmin_meanAnnAvg_CLIM,
"tmax" = tmax_meanAnnAvg_CLIM, #1
"tmean" = tmean_meanAnnAvg_CLIM,
"prcp" = prcp_meanAnnTotal_CLIM,
"t_warm" = T_warmestMonth_meanAnnAvg_CLIM,
"t_cold" = T_coldestMonth_meanAnnAvg_CLIM,
"prcp_wet" = precip_wettestMonth_meanAnnAvg_CLIM,
"prcp_dry" = precip_driestMonth_meanAnnAvg_CLIM,
"prcp_seasonality" = precip_Seasonality_meanAnnAvg_CLIM, #2
"prcpTempCorr" = PrecipTempCorr_meanAnnAvg_CLIM, #3
"abvFreezingMonth" = aboveFreezing_month_meanAnnAvg_CLIM,
"isothermality" = isothermality_meanAnnAvg_CLIM, #4
"annWatDef" = annWaterDeficit_meanAnnAvg_CLIM,
"annWetDegDays" = annWetDegDays_meanAnnAvg_CLIM,
"VPD_mean" = annVPD_mean_meanAnnAvg_CLIM,
"VPD_max" = annVPD_max_meanAnnAvg_CLIM, #5
"VPD_min" = annVPD_min_meanAnnAvg_CLIM, #6
"VPD_max_95" = annVPD_max_95percentile_CLIM,
"annWatDef_95" = annWaterDeficit_95percentile_CLIM,
"annWetDegDays_5" = annWetDegDays_5percentile_CLIM,
"frostFreeDays_5" = durationFrostFreeDays_5percentile_CLIM,
"frostFreeDays" = durationFrostFreeDays_meanAnnAvg_CLIM,
"soilDepth" = soilDepth, #7
"clay" = surfaceClay_perc,
"sand" = avgSandPerc_acrossDepth, #8
"coarse" = avgCoarsePerc_acrossDepth, #9
"carbon" = avgOrganicCarbonPerc_0_3cm, #10
"AWHC" = totalAvailableWaterHoldingCapacity,
## anomaly variables
tmean_anom = tmean_meanAnnAvg_3yrAnom, #15
tmin_anom = tmin_meanAnnAvg_3yrAnom, #16
tmax_anom = tmax_meanAnnAvg_3yrAnom, #17
prcp_anom = prcp_meanAnnTotal_3yrAnom, #18
t_warm_anom = T_warmestMonth_meanAnnAvg_3yrAnom, #19
t_cold_anom = T_coldestMonth_meanAnnAvg_3yrAnom, #20
prcp_wet_anom = precip_wettestMonth_meanAnnAvg_3yrAnom, #21
precp_dry_anom = precip_driestMonth_meanAnnAvg_3yrAnom, #22
prcp_seasonality_anom = precip_Seasonality_meanAnnAvg_3yrAnom, #23
prcpTempCorr_anom = PrecipTempCorr_meanAnnAvg_3yrAnom, #24
aboveFreezingMonth_anom = aboveFreezing_month_meanAnnAvg_3yrAnom, #25
isothermality_anom = isothermality_meanAnnAvg_3yrAnom, #26
annWatDef_anom = annWaterDeficit_meanAnnAvg_3yrAnom, #27
annWetDegDays_anom = annWetDegDays_meanAnnAvg_3yrAnom, #28
VPD_mean_anom = annVPD_mean_meanAnnAvg_3yrAnom, #29
VPD_min_anom = annVPD_min_meanAnnAvg_3yrAnom, #30
VPD_max_anom = annVPD_max_meanAnnAvg_3yrAnom, #31
VPD_max_95_anom = annVPD_max_95percentile_3yrAnom, #32
annWatDef_95_anom = annWaterDeficit_95percentile_3yrAnom, #33
annWetDegDays_5_anom = annWetDegDays_5percentile_3yrAnom , #34
frostFreeDays_5_anom = durationFrostFreeDays_5percentile_3yrAnom, #35
frostFreeDays_anom = durationFrostFreeDays_meanAnnAvg_3yrAnom #36
)
# small dataset for if testing the data
if(test_run) {
modDat_1 <- slice_sample(modDat_1, n = 1e5)
}Add a constant to the response variable (+1) so that models run…
modDat_1[,response] <- modDat_1[,response]+2Identify the ecoregion and response variable type to use in this model run
ecoregion <- params$ecoregion
response <- params$response
print(paste0("In this model run, the ecoregion is ", ecoregion," and the response variable is ",response))## [1] "In this model run, the ecoregion is shrubGrass and the response variable is BareGroundCover"Subset the data to only include data for the ecoregion of interest
if (ecoregion == "shrubGrass") {
# select data for the ecoregion of interest
modDat_1 <- modDat_1 %>%
filter(newRegion == "dryShrubGrass")
} else if (ecoregion == "forest") {
# select data for the ecoregion of interest
modDat_1 <- modDat_1 %>%
filter(newRegion %in% c("eastForest", "westForest"))
}
# remove the rows that have no observations for the response variable of interest
modDat_1 <- modDat_1[!is.na(modDat_1[,response]),]Currently, subsampling data from the “Texas Coastal Plain”, since it’s quite different from other regions and is really messing with model fit
modDat_1_noLA <- modDat_1 %>%
filter(NA_L2NAME != "TEXAS-LOUISIANA COASTAL PLAIN")
modDat_1_LA <- modDat_1 %>%
filter(NA_L2NAME == "TEXAS-LOUISIANA COASTAL PLAIN")
# sample points
modDat_1 <- modDat_1_LA %>%
slice_sample(n = round(nrow(modDat_1_LA)*.3)) %>%
rbind(modDat_1_noLA) Visualize the response variable
hist(modDat_1[,response], main = paste0("Histogram of ",response),
xlab = paste0(response))
Visualize the predictor variables
The following are the candidate predictor variables for this ecoregion:
if (ecoregion == "shrubGrass") {
# select potential predictor variables for the ecoregion of interest
prednames <-
c(
"tmean" , "prcp" ,"prcp_seasonality" ,"prcpTempCorr" ,
"isothermality" , "annWatDef" ,"sand" ,"coarse" ,
"carbon" , "AWHC" ,"tmin_anom" ,"tmax_anom" ,
"t_warm_anom" , "prcp_wet_anom" ,"precp_dry_anom" ,"prcp_seasonality_anom" ,
"prcpTempCorr_anom" , "aboveFreezingMonth_anom" ,"isothermality_anom" ,"annWatDef_anom" ,
"annWetDegDays_anom", "VPD_mean_anom" ,"VPD_min_anom" ,"frostFreeDays_5_anom" )
} else if (ecoregion == "forest") {
# select potential predictor variables for the ecoregion of interest
prednames <-
c(
"tmean" ,"prcp" , "prcp_dry" , "prcpTempCorr" ,
"isothermality" ,"annWatDef" , "clay" , "sand" ,
"coarse" ,"carbon" , "AWHC" , "tmin_anom" ,
"tmax_anom" ,"prcp_anom" , "prcp_wet_anom" , "precp_dry_anom" ,
"prcp_seasonality_anom" ,"prcpTempCorr_anom" , "aboveFreezingMonth_anom" , "isothermality_anom",
"annWatDef_anom" ,"annWetDegDays_anom" , "VPD_mean_anom" , "VPD_max_95_anom" ,
"frostFreeDays_5_anom" )
}
# subset the data to only include these predictors, and remove any remaining NAs
modDat_1 <- modDat_1 %>%
select(prednames, response, newRegion, Year, Long, Lat, NA_L1NAME, NA_L2NAME) %>%
drop_na()
names(prednames) <- prednames
df_pred <- modDat_1[, prednames]
#
# # print the list of predictor variables
# knitr::kable(format = "html", data.frame("Possible_Predictors" = prednames)
# ) %>%
# kable_styling(bootstrap_options = c("striped", "hover", "condensed")) create_summary <- function(df) {
df %>%
pivot_longer(cols = everything(),
names_to = 'variable') %>%
group_by(variable) %>%
summarise(across(value, .fns = list(mean = ~mean(.x, na.rm = TRUE), min = ~min(.x, na.rm = TRUE),
median = ~median(.x, na.rm = TRUE), max = ~max(.x, na.rm = TRUE)))) %>%
mutate(across(where(is.numeric), round, 4))
}
modDat_1[prednames] %>%
create_summary() %>%
knitr::kable(caption = 'summaries of possible predictor variables') %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed")) | variable | value_mean | value_min | value_median | value_max |
|---|---|---|---|---|
| AWHC | 14.4682 | 0.6319 | 13.9491 | 35.1058 |
| VPD_mean_anom | -0.0191 | -0.3533 | -0.0181 | 0.2146 |
| VPD_min_anom | -0.0118 | -0.1747 | -0.0147 | 0.1390 |
| aboveFreezingMonth_anom | 0.0576 | -1.3667 | 0.0417 | 2.2143 |
| annWatDef | 99.6331 | 1.5871 | 88.3355 | 451.1875 |
| annWatDef_anom | -0.0591 | -3.9347 | -0.0570 | 1.0000 |
| annWetDegDays_anom | 0.0023 | -1.6606 | 0.0037 | 0.9840 |
| carbon | 1.0402 | 0.0001 | 0.8316 | 27.0162 |
| coarse | 8.3024 | 0.0000 | 5.4085 | 63.1943 |
| frostFreeDays_5_anom | -16.6301 | -216.0000 | -6.0000 | 55.8000 |
| isothermality | 37.4422 | 23.2374 | 37.2059 | 54.8325 |
| isothermality_anom | 0.4858 | -6.0230 | 0.4179 | 8.9119 |
| prcp | 417.0681 | 90.8317 | 371.8310 | 1695.6347 |
| prcpTempCorr | 0.0531 | -0.7694 | 0.1101 | 0.6890 |
| prcpTempCorr_anom | 0.0170 | -0.6495 | 0.0194 | 0.5744 |
| prcp_seasonality | 0.9366 | 0.5255 | 0.9216 | 1.7468 |
| prcp_seasonality_anom | -0.0138 | -0.8313 | -0.0061 | 0.4743 |
| prcp_wet_anom | -0.0131 | -1.9576 | 0.0068 | 0.6958 |
| precp_dry_anom | -0.0547 | -9.0000 | 0.2449 | 1.0000 |
| sand | 45.9082 | 2.5763 | 45.5166 | 98.9991 |
| t_warm_anom | -0.4440 | -6.3974 | -0.4584 | 3.3598 |
| tmax_anom | -0.2635 | -3.4286 | -0.2991 | 2.3198 |
| tmean | 10.5113 | -2.2935 | 9.2919 | 23.9289 |
| tmin_anom | -0.4478 | -4.0550 | -0.4295 | 2.4361 |
# response_summary <- modDat_1 %>%
# dplyr::select(#where(is.numeric), -all_of(pred_vars),
# matches(response)) %>%
# create_summary()
#
#
# kable(response_summary,
# caption = 'summaries of response variables, calculated using paint') %>%
# kable_styling(bootstrap_options = c("striped", "hover", "condensed")) Plot predictor vars against each other
set.seed(12011993)
# function for colors
my_fn <- function(data, mapping, method="p", use="pairwise", ...){
# grab data
x <- eval_data_col(data, mapping$x)
y <- eval_data_col(data, mapping$y)
# calculate correlation
corr <- cor(x, y, method=method, use=use)
# calculate colour based on correlation value
# Here I have set a correlation of minus one to blue,
# zero to white, and one to red
# Change this to suit: possibly extend to add as an argument of `my_fn`
colFn <- colorRampPalette(c("red", "white", "blue"), interpolate ='spline')
fill <- colFn(100)[findInterval(corr, seq(-1, 1, length=100))]
ggally_cor(data = data, mapping = mapping, size = 2.5, stars = FALSE,
digits = 2, colour = I("black"),...) +
theme_void() +
theme(panel.background = element_rect(fill=fill))
}
if (run == TRUE) {
(corrPlot <- modDat_1 %>%
select(prednames) %>%
slice_sample(n = 5e4) %>%
#select(-matches("_")) %>%
ggpairs( upper = list(continuous = my_fn, size = .1), lower = list(continuous = GGally::wrap("points", alpha = 0.1, size=0.1)), progress = FALSE))
base::saveRDS(corrPlot, paste0("../ModelFitting/models/", response, "_",ecoregion, "_corrPlot.rds"))
} else {
# corrPlot <- readRDS(paste0("../ModelFitting/models/", response, "_",ecoregion, "_corrPlot.rds"))
# (corrPlot)
print(c("See previous correlation figures"))
}Predictor variables compared to binned response variables
set.seed(12011993)
# vector of name of response variables
vars_response <- response
# longformat dataframes for making boxplots
df_sample_plots <- modDat_1 %>%
slice_sample(n = 5e4) %>%
rename(response = all_of(response)) %>%
mutate(response = case_when(
response <= .25 ~ ".25",
response > .25 & response <=.5 ~ ".5",
response > .5 & response <=.75 ~ ".75",
response >= .75 ~ "1",
)) %>%
select(c(response, prednames)) %>%
tidyr::pivot_longer(cols = unname(prednames),
names_to = "predictor",
values_to = "value"
)
ggplot(df_sample_plots, aes_string(x= "response", y = 'value')) +
geom_boxplot() +
facet_wrap(~predictor , scales = 'free_y') +
ylab("Predictor Variable Values") +
xlab(response)Standardize the predictor variables for the model-fitting process
modDat_1_s <- modDat_1 %>%
mutate(across(all_of(prednames), base::scale, .names = "{.col}_s"))
names(modDat_1_s) <- c(names(modDat_1),
paste0(prednames, "_s")
)
scaleFigDat_1 <- modDat_1_s %>%
select(c(Long, Lat, Year, prednames)) %>%
pivot_longer(cols = all_of(names(prednames)),
names_to = "predNames",
values_to = "predValues_unScaled")
scaleFigDat_2 <- modDat_1_s %>%
select(c(Long, Lat, Year,paste0(prednames, "_s"))) %>%
pivot_longer(cols = all_of(paste0(prednames, "_s")),
names_to = "predNames",
values_to = "predValues_scaled",
names_sep = ) %>%
mutate(predNames = str_replace(predNames, pattern = "_s$", replacement = ""))
scaleFigDat_3 <- scaleFigDat_1 %>%
left_join(scaleFigDat_2)
ggplot(scaleFigDat_3) +
facet_wrap(~predNames, scales = "free") +
geom_histogram(aes(predValues_unScaled), fill = "lightgrey", col = "darkgrey") +
geom_histogram(aes(predValues_scaled), fill = "lightblue", col = "blue") +
xlab ("predictor variable values") +
ggtitle("Comparing the distribution of unscaled (grey) to scaled (blue) predictor variables")
Model Fitting
Visualize the level 2 ecoregions and how they differ across environmental space
## visualize the variation between groups across environmental space
## do a pca of climate across level 2 ecoregions
pca <- prcomp(modDat_1_s[,paste0(prednames, "_s")])
library(factoextra)
(fviz_pca_ind(pca, habillage = modDat_1_s$NA_L2NAME, label = "none", addEllipses = TRUE, ellipse.level = .95, ggtheme = theme_minimal(), alpha.ind = .1))
if (ecoregion == "shrubGrass") {
print("We'll combine the 'Mediterranean California' and 'Western Sierra Madre Piedmont' ecoregions (into 'Mediterranean Piedmont'). We'll also combine `Tamaulipas-Texas semiarid plain,' 'Texas-Lousiana Coastal plain,' and 'South Central semiarid prairies' ecoregions (into (`Semiarid plain and prairies`)." )
modDat_1_s[modDat_1_s$NA_L2NAME %in% c("MEDITERRANEAN CALIFORNIA", "WESTERN SIERRA MADRE PIEDMONT"), "NA_L2NAME"] <- "MEDITERRANEAN PIEDMONT"
modDat_1[modDat_1$NA_L2NAME %in% c("MEDITERRANEAN CALIFORNIA", "WESTERN SIERRA MADRE PIEDMONT"), "NA_L2NAME"] <- "MEDITERRANEAN PIEDMONT"
modDat_1_s[modDat_1_s$NA_L2NAME %in% c("TAMAULIPAS-TEXAS SEMIARID PLAIN", "TEXAS-LOUISIANA COASTAL PLAIN", "SOUTH CENTRAL SEMIARID PRAIRIES"), "NA_L2NAME"] <- "SEMIARID PLAIN AND PRAIRIES"
modDat_1[modDat_1$NA_L2NAME %in% c("TAMAULIPAS-TEXAS SEMIARID PLAIN", "TEXAS-LOUISIANA COASTAL PLAIN", "SOUTH CENTRAL SEMIARID PRAIRIES"), "NA_L2NAME"] <- "SEMIARID PLAIN AND PRAIRIES"
}## [1] "We'll combine the 'Mediterranean California' and 'Western Sierra Madre Piedmont' ecoregions (into 'Mediterranean Piedmont'). We'll also combine `Tamaulipas-Texas semiarid plain,' 'Texas-Lousiana Coastal plain,' and 'South Central semiarid prairies' ecoregions (into (`Semiarid plain and prairies`)."# make a table of n for each region
modDat_1 %>%
group_by(NA_L2NAME) %>%
dplyr::summarize("Number_Of_Observations" = length(NA_L2NAME)) %>%
rename("Level_2_Ecoregion" = NA_L2NAME)%>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed")) | Level_2_Ecoregion | Number_Of_Observations |
|---|---|
| COLD DESERTS | 63161 |
| MEDITERRANEAN PIEDMONT | 2606 |
| SEMIARID PLAIN AND PRAIRIES | 24921 |
| TEMPERATE PRAIRIES | 3045 |
| WARM DESERTS | 10353 |
| WEST-CENTRAL SEMIARID PRAIRIES | 19493 |
Then, look at the spatial distribution and environmental characteristics of the grouped ecoregions
## make data into spatial format
modDat_1_sf <- modDat_1 %>%
st_as_sf(coords = c("Long", "Lat"), crs = st_crs("PROJCRS[\"unnamed\",\n BASEGEOGCRS[\"unknown\",\n DATUM[\"unknown\",\n ELLIPSOID[\"Spheroid\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433,\n ID[\"EPSG\",9122]]]],\n CONVERSION[\"Lambert Conic Conformal (2SP)\",\n METHOD[\"Lambert Conic Conformal (2SP)\",\n ID[\"EPSG\",9802]],\n PARAMETER[\"Latitude of false origin\",42.5,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8821]],\n PARAMETER[\"Longitude of false origin\",-100,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8822]],\n PARAMETER[\"Latitude of 1st standard parallel\",25,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8823]],\n PARAMETER[\"Latitude of 2nd standard parallel\",60,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8824]],\n PARAMETER[\"Easting at false origin\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8826]],\n PARAMETER[\"Northing at false origin\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8827]]],\n CS[Cartesian,2],\n AXIS[\"easting\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]],\n AXIS[\"northing\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]]]"))
# download map info for visualization
data(state_boundaries_wgs84)
cropped_states <- suppressMessages(state_boundaries_wgs84 %>%
dplyr::filter(NAME!="Hawaii") %>%
dplyr::filter(NAME!="Alaska") %>%
dplyr::filter(NAME!="Puerto Rico") %>%
dplyr::filter(NAME!="American Samoa") %>%
dplyr::filter(NAME!="Guam") %>%
dplyr::filter(NAME!="Commonwealth of the Northern Mariana Islands") %>%
dplyr::filter(NAME!="United States Virgin Islands") %>%
sf::st_sf() %>%
sf::st_transform(sf::st_crs(modDat_1_sf))) #%>%
#sf::st_crop(sf::st_bbox(modDat_1_sf)+c(-1,-1,1,1))
map1 <- ggplot() +
geom_sf(data=cropped_states,fill='white') +
geom_sf(data=modDat_1_sf,aes(fill=as.factor(NA_L2NAME)),linewidth=0.5,alpha=0.5) +
geom_point(data=modDat_1,alpha=0.5,
aes(x = Long, y = Lat, color=as.factor(NA_L2NAME)), alpha = .1) +
#scale_fill_okabeito() +
#scale_color_okabeito() +
# theme_default() +
theme(legend.position = 'none') +
labs(title = "Level 2 Ecoregions as spatial blocks")
hull <- modDat_1_sf %>%
ungroup() %>%
group_by(NA_L2NAME) %>%
slice(chull(tmean, prcp))
plot1<-ggplot(data=modDat_1_sf,aes(x=tmean,y=prcp)) +
geom_polygon(data = hull, alpha = 0.25,aes(fill=NA_L2NAME) )+
geom_point(aes(group=NA_L2NAME,color=NA_L2NAME),alpha=0.25) +
theme_minimal() + xlab("Annual Average T_mean - long-term average") +
ylab("Annual Average Precip - long-term average") #+
#scale_color_okabeito() +
#scale_fill_okabeito()
plot2<-ggplot(data=modDat_1_sf %>%
pivot_longer(cols=tmean:prcp),
aes(x=value,group=name)) +
# geom_polygon(data = hull, alpha = 0.25,aes(fill=fold) )+
geom_density(aes(group=NA_L2NAME,fill=NA_L2NAME),alpha=0.25) +
theme_minimal() +
facet_wrap(~name,scales='free')# +
#scale_color_okabeito() +
#scale_fill_okabeito()
library(patchwork)
(combo <- (map1+plot1)/plot2) 
Fit a global model with all of the data
First, fit a LASSO regression model using the glmnet R package
- This regression is a Gamma glm with a log link
- Use cross validation across level 2 ecoregions to tune the lambda parameter in the LASSO model
- this model is fit to using the scaled weather/climate/soils variables
- this list of possible predictors includes:
- main effects
- interactions between all soils variables
- interactions between climate and weather variables
- transformed main effects (squared, log-transformed (add a uniform integer – 20– to all variables prior to log-transformation), square root-transformed (add a uniform integer – 20– to all variables prior to log-transformation))
##
## Call: cv.glmnet(x = X[, 2:ncol(X)], y = y, type.measure = "mse", foldid = my_folds, keep = TRUE, parallel = TRUE, family = stats::Gamma(link = "log"), alpha = 1, nlambda = 100, standardize = FALSE)
##
## Measure: Mean-Squared Error
##
## Lambda Index Measure SE Nonzero
## min 0.00579 49 210.5 43.60 78
## 1se 0.05402 25 249.8 53.69 15
Then, fit regular glm models (Gamma glm with a log link), first using the coefficients from the ‘best’ lambda identified in the LASSO model, as then using the coefficients from the ‘1SE’ lambda identified from the LASSO (this is the value of lambda such that the cross validation error is within 1 standard error of the minimum).
## fit w/ the identified coefficients from the 'best' lambda, but using the glm function
mat_glmnet_best <- as.matrix(bestLambda_coef)
mat2_glmnet_best <- mat_glmnet_best[mat_glmnet_best[,1] != 0,]
names(mat2_glmnet_best) <- rownames(mat_glmnet_best)[mat_glmnet_best[,1] != 0]
if (length(mat2_glmnet_best) == 1) {
f_glm_bestLambda <- as.formula(paste0(response, "~ 1"))
} else {
f_glm_bestLambda <- as.formula(paste0(response, " ~ ", paste0(names( mat2_glmnet_best)[2:length(names( mat2_glmnet_best))], collapse = " + ")))
}
fit_glm_bestLambda <- glm(data = modDat_1_s
, formula = f_glm_bestLambda, family = stats::Gamma(link = "log"))
## fit w/ the identified coefficients from the '1se' lambda, but using the glm function
mat_glmnet_1se <- as.matrix(seLambda_coef)
mat2_glmnet_1se <- mat_glmnet_1se[mat_glmnet_1se[,1] != 0,]
names(mat2_glmnet_1se) <- rownames(mat_glmnet_1se)[mat_glmnet_1se[,1] != 0]
if(length(mat2_glmnet_1se) == 1) {
f_glm_1se <- as.formula(paste0(response, "~ 1"))
} else {
f_glm_1se <- as.formula(paste0(response, " ~ ", paste0(names( mat2_glmnet_1se)[2:length(names( mat2_glmnet_1se))], collapse = " + ")))
}
fit_glm_se <- glm(data = modDat_1_s, formula = f_glm_1se,
family = stats::Gamma(link = "log"))Then, we predict (on the training set) using both of these models (best lambda and 1se lambda)
## predict on the test data
# lasso model predictions with the optimal lambda
optimal_pred <- predict(fit_glm_bestLambda, newx=X[,2:ncol(X)], type = "response") - 2
optimal_pred_1se <- predict(fit_glm_se, newx=X[,2:ncol(X)], type = "response") - 2
null_fit <- glm(#data = data.frame("y" = y, X[,paste0(prednames, "_s")]),
formula = y ~ 1, family = stats::Gamma(link = "log"))
null_pred <- predict(null_fit, newdata = as.data.frame(X), type = "response"
) - 2
# save data
fullModOut <- list(
"modelObject" = fit,
"nullModelObject" = null_fit,
"modelPredictions" = data.frame(#ecoRegion_holdout = rep(test_eco,length(y)),
obs=y,
pred_opt=optimal_pred,
pred_opt_se = optimal_pred_1se,
pred_null=null_pred#,
#pred_nopenalty=nopen_pred
))
# calculate correlations between null and optimal model
my_cors <- c(cor(optimal_pred, c(y-2)),
cor(optimal_pred_1se, c(y-2)),
cor(null_pred, c(y-2))
)
# calculate mse between null and optimal model
my_mse <- c(mean((fullModOut$modelPredictions$pred_opt - c(y-2))^2) ,
mean((fullModOut$modelPredictions$pred_opt_se - c(y-2))^2) ,
mean((fullModOut$modelPredictions$pred_null - c(y-2))^2)#,
#mean((obs_pred$pred_nopenalty - obs_pred$obs)^2)
)
ggplot() +
geom_point(aes(X[,2], fullModOut$modelPredictions$obs-2), col = "black", alpha = .1) +
geom_point(aes(X[,2], fullModOut$modelPredictions$pred_opt), col = "red", alpha = .1) + ## predictions w/ the CV model
geom_point(aes(X[,2], fullModOut$modelPredictions$pred_opt_se), col = "green", alpha = .1) + ## predictions w/ the CV model (1se lambda)
geom_point(aes(X[,2], fullModOut$modelPredictions$pred_null), col = "blue", alpha = .1) +
labs(title = "A rough comparison of observed and model-predicted values",
subtitle = "black = observed values \n red = predictions from 'best lambda' model \n green = predictions from '1se' lambda model \n blue = predictions from null model") +
xlab(colnames(X)[2])
#ylim(c(0,200))The internal cross-validation process to fit the global LASSO model
identified an optimal lambda value (regularization parameter) of
r{print(best_lambda)}. The lambda value such that the cross
validation error is within 1 standard error of the minimum (“1se
lambda”) was `r{print(fit$lambda.1se)}`` . The following coefficients
were kept in each model:
# the coefficient matrix from the 'best model' -- find and print those coefficients that aren't 0 in a table
coef_glm_bestLambda <- coef(fit_glm_bestLambda) %>%
data.frame()
coef_glm_bestLambda$coefficientName <- rownames(coef_glm_bestLambda)
names(coef_glm_bestLambda)[1] <- "coefficientValue_bestLambda"
# coefficient matrix from the '1se' model
coef_glm_1se <- coef(fit_glm_se) %>%
data.frame()
coef_glm_1se$coefficientName <- rownames(coef_glm_1se)
names(coef_glm_1se)[1] <- "coefficientValue_1seLambda"
# add together
coefs <- full_join(coef_glm_bestLambda, coef_glm_1se) %>%
select(coefficientName, coefficientValue_bestLambda, coefficientValue_1seLambda)
globModTerms <- coefs[!is.na(coefs$coefficientValue_bestLambda), "coefficientName"]
## also, get the number of unique variables in each model
var_prop_pred <- paste0(response, "_pred")
response_vars <- c(response, var_prop_pred)
# for best lambda model
prednames_fig <- paste(str_split(globModTerms, ":", simplify = TRUE))
prednames_fig <- str_replace(prednames_fig, "I\\(", "")
prednames_fig <- str_replace(prednames_fig, "\\^2\\)", "")
prednames_fig <- unique(prednames_fig[prednames_fig>0])
prednames_fig <- prednames_fig
prednames_fig_num <- length(prednames_fig)
# for 1SE lambda model
globModTerms_1se <- coefs[!is.na(coefs$coefficientValue_1seLambda), "coefficientName"]
if (length(globModTerms_1se) == 1) {
prednames_fig_1se <- paste(str_split(globModTerms_1se, ":", simplify = TRUE))
prednames_fig_1se <- str_replace(prednames_fig_1se, "I\\(", "")
prednames_fig_1se <- str_replace(prednames_fig_1se, "\\^2\\)", "")
prednames_fig_1se <- unique(prednames_fig_1se[prednames_fig_1se>0])
prednames_fig_1se_num <- c(0)
} else {
prednames_fig_1se <- paste(str_split(globModTerms_1se, ":", simplify = TRUE))
prednames_fig_1se <- str_replace(prednames_fig_1se, "I\\(", "")
prednames_fig_1se <- str_replace(prednames_fig_1se, "\\^2\\)", "")
prednames_fig_1se <- unique(prednames_fig_1se[prednames_fig_1se>0])
prednames_fig_1se_num <- length(prednames_fig_1se)
}
# make a table
kable(coefs, col.names = c("Coefficient Name", "Value from best lambda model", "Value from 1se lambda model")
) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed")) | Coefficient Name | Value from best lambda model | Value from 1se lambda model |
|---|---|---|
| (Intercept) | 2.7607260 | 2.8133793 |
| prcp_s | -0.6486547 | -0.5755896 |
| prcp_seasonality_s | -0.1008334 | NA |
| prcpTempCorr_s | -0.0594209 | -0.0857223 |
| isothermality_s | 0.1437410 | 0.1277627 |
| sand_s | -0.0374452 | NA |
| coarse_s | -0.0800781 | NA |
| carbon_s | -0.2656310 | -0.2692416 |
| t_warm_anom_s | -0.0350100 | -0.0351340 |
| prcp_wet_anom_s | 0.0525289 | NA |
| prcp_seasonality_anom_s | -0.0420171 | NA |
| aboveFreezingMonth_anom_s | -0.0023656 | NA |
| isothermality_anom_s | -0.0077912 | NA |
| annWatDef_anom_s | -0.0007452 | -0.0332761 |
| VPD_mean_anom_s | -0.0199192 | NA |
| I(prcp_s^2) | 0.0975047 | 0.0760063 |
| I(annWatDef_s^2) | -0.0184579 | NA |
| I(t_warm_anom_s^2) | -0.0089036 | NA |
| I(prcp_wet_anom_s^2) | 0.0022478 | NA |
| I(precp_dry_anom_s^2) | -0.0063902 | -0.0051332 |
| I(aboveFreezingMonth_anom_s^2) | -0.0033975 | NA |
| I(isothermality_anom_s^2) | -0.0038881 | NA |
| I(annWatDef_anom_s^2) | 0.0057781 | 0.0036480 |
| I(VPD_mean_anom_s^2) | 0.0146675 | NA |
| I(sand_s^2) | 0.0260941 | NA |
| I(carbon_s^2) | 0.0162096 | 0.0133649 |
| I(AWHC_s^2) | 0.0049890 | NA |
| prcp_wet_anom_s:aboveFreezingMonth_anom_s | 0.0072377 | NA |
| aboveFreezingMonth_anom_s:prcpTempCorr_anom_s | -0.0056280 | NA |
| isothermality_s:annWatDef_anom_s | -0.0093516 | NA |
| annWatDef_anom_s:tmean_s | -0.0205678 | NA |
| annWatDef_anom_s:VPD_mean_anom_s | -0.0137775 | NA |
| isothermality_s:annWatDef_s | -0.0378093 | NA |
| prcpTempCorr_s:annWatDef_s | 0.0676267 | 0.1274036 |
| annWatDef_s:precp_dry_anom_s | -0.0156467 | NA |
| isothermality_anom_s:annWetDegDays_anom_s | -0.0036792 | NA |
| prcp_wet_anom_s:annWetDegDays_anom_s | 0.0110703 | NA |
| isothermality_s:frostFreeDays_5_anom_s | 0.0161090 | NA |
| prcpTempCorr_s:frostFreeDays_5_anom_s | -0.0133376 | NA |
| t_warm_anom_s:frostFreeDays_5_anom_s | 0.0019398 | NA |
| frostFreeDays_5_anom_s:tmax_anom_s | -0.0078261 | NA |
| isothermality_s:isothermality_anom_s | -0.0179137 | NA |
| prcp_s:isothermality_anom_s | -0.0203626 | NA |
| prcp_seasonality_s:isothermality_anom_s | 0.0106286 | NA |
| prcpTempCorr_s:isothermality_anom_s | 0.0024360 | NA |
| t_warm_anom_s:isothermality_anom_s | 0.0111584 | NA |
| prcp_s:isothermality_s | -0.0454616 | NA |
| prcp_seasonality_s:isothermality_s | -0.0270188 | NA |
| isothermality_s:prcp_wet_anom_s | -0.0093031 | NA |
| prcpTempCorr_s:isothermality_s | -0.0457330 | -0.0022826 |
| isothermality_s:tmax_anom_s | -0.0142512 | NA |
| prcp_s:prcp_wet_anom_s | 0.0092251 | NA |
| prcp_s:prcpTempCorr_s | -0.1200517 | NA |
| prcp_s:t_warm_anom_s | 0.0075792 | NA |
| prcp_s:tmean_s | -0.0503040 | NA |
| prcp_s:VPD_mean_anom_s | 0.0121169 | NA |
| prcp_seasonality_anom_s:prcpTempCorr_anom_s | 0.0015387 | NA |
| prcpTempCorr_s:prcp_seasonality_anom_s | -0.0062247 | NA |
| t_warm_anom_s:prcp_seasonality_anom_s | -0.0159484 | NA |
| prcp_seasonality_anom_s:VPD_mean_anom_s | 0.0132898 | NA |
| prcp_seasonality_s:prcp_wet_anom_s | -0.0049159 | NA |
| prcp_seasonality_s:prcpTempCorr_s | -0.0653724 | NA |
| prcp_seasonality_s:tmean_s | 0.0570051 | NA |
| prcp_seasonality_s:VPD_min_anom_s | 0.0074978 | NA |
| prcp_wet_anom_s:tmax_anom_s | -0.0051103 | NA |
| prcpTempCorr_s:prcpTempCorr_anom_s | 0.0197944 | NA |
| prcpTempCorr_s:tmean_s | 0.1085326 | NA |
| prcpTempCorr_s:tmin_anom_s | -0.0201332 | NA |
| prcpTempCorr_s:VPD_mean_anom_s | 0.0154403 | NA |
| t_warm_anom_s:VPD_mean_anom_s | -0.0126546 | NA |
| tmean_s:tmax_anom_s | -0.0100774 | NA |
| VPD_mean_anom_s:tmean_s | 0.0072887 | NA |
| tmean_s:VPD_min_anom_s | 0.0080997 | NA |
| VPD_mean_anom_s:VPD_min_anom_s | -0.0218783 | NA |
| carbon_s:AWHC_s | -0.0588676 | NA |
| coarse_s:AWHC_s | 0.0120734 | NA |
| sand_s:AWHC_s | 0.0198336 | NA |
| coarse_s:carbon_s | 0.0207331 | NA |
| sand_s:carbon_s | 0.0113956 | NA |
| I(isothermality_s^2) | NA | -0.0182360 |
| I(coarse_s^2) | NA | -0.0255511 |
| carbon_s:coarse_s | NA | 0.0625665 |
# calculate RMSE of both models
RMSE_best <- yardstick::rmse(fullModOut$modelPredictions[,c("obs", "pred_opt")]-2, truth = "obs", estimate = "pred_opt")$.estimate
RMSE_1se <- yardstick::rmse(fullModOut$modelPredictions[,c("obs", "pred_opt_se")]-2, truth = "obs", estimate = "pred_opt_se")$.estimate
bias_best <- mean(fullModOut$modelPredictions$obs - fullModOut$modelPredictions$pred_opt)
bias_1se <- mean(fullModOut$modelPredictions$obs - fullModOut$modelPredictions$pred_opt_se)
uniqueCoeffs <- data.frame("Best lambda model" = c(RMSE_best, bias_best,
as.integer(length(globModTerms)-1), as.integer(prednames_fig_num),
as.integer(sum(prednames_fig %in% c(prednames_clim))),
as.integer(sum(prednames_fig %in% c(prednames_weath))),
as.integer(sum(prednames_fig %in% c(prednames_soils)))
),
"1se lambda model" = c(RMSE_1se, bias_1se,
length(globModTerms_1se)-1, prednames_fig_1se_num,
sum(prednames_fig_1se %in% c(prednames_clim)),
sum(prednames_fig_1se %in% c(prednames_weath)),
sum(prednames_fig_1se %in% c(prednames_soils))))
row.names(uniqueCoeffs) <- c("RMSE", "bias - mean(obs-pred.)", "Total number of coefficients", "Number of unique coefficients",
"Number of unique climate coefficients",
"Number of unique weather coefficients",
"Number of unique soils coefficients"
)
kable(uniqueCoeffs,
col.names = c("Best lambda model", "1se lambda model"), row.names = TRUE) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed")) | Best lambda model | 1se lambda model | |
|---|---|---|
| RMSE | 12.978342 | 13.488516 |
| bias - mean(obs-pred.) | 2.132895 | 2.138661 |
| Total number of coefficients | 78.000000 | 15.000000 |
| Number of unique coefficients | 24.000000 | 9.000000 |
| Number of unique climate coefficients | 6.000000 | 4.000000 |
| Number of unique weather coefficients | 14.000000 | 3.000000 |
| Number of unique soils coefficients | 4.000000 | 2.000000 |
Visualizations of Model Predictions and Residuals – using best lambda model
observed vs. predicted values
Predicting on the data
# create prediction for each each model
# (i.e. for each fire proporation variable)
predict_by_response <- function(mod, df) {
df_out <- df
response_name <- paste0(response, "_pred")
df_out <- df_out %>% cbind(predict(mod, newx= df_out, #s="lambda.min",
type = "response"))
colnames(df_out)[ncol(df_out)] <- response_name
return(df_out)
}
pred_glm1 <- predict_by_response(fit_glm_bestLambda, X[,2:ncol(X)])
# add back in true y values
pred_glm1 <- pred_glm1 %>%
cbind( data.frame("y" = y))
# rename the true response column to not be 'y_test'
colnames(pred_glm1)[which(colnames(pred_glm1) == "y")] <- paste(response)
# add back in lat/long data
pred_glm1 <- pred_glm1 %>%
cbind(modDat_1_s[,c("Long", "Lat", "Year")])
pred_glm1$resid <- pred_glm1[,response] - pred_glm1[,paste0(response, "_pred")]
pred_glm1$extremeResid <- NA
pred_glm1[pred_glm1$resid > 70 | pred_glm1$resid < -70,"extremeResid"] <- 1
# plot(x = pred_glm1[,response],
# y = pred_glm1[,paste0(response, "_pred")],
# xlab = "observed values", ylab = "predicted values")
# points(x = pred_glm1[!is.na(pred_glm1$extremeResid), response],
# y = pred_glm1[!is.na(pred_glm1$extremeResid), paste0(response, "_pred")],
# col = "red")pred_glm1_1se <- predict_by_response(fit_glm_se, X[,2:ncol(X)])
# add back in true y values
pred_glm1_1se <- pred_glm1_1se %>%
cbind( data.frame("y" = y))
# rename the true response column to not be 'y_test'
colnames(pred_glm1_1se)[which(colnames(pred_glm1_1se) == "y")] <- paste(response)
# add back in lat/long data
pred_glm1_1se <- pred_glm1_1se %>%
cbind(modDat_1_s[,c("Long", "Lat", "Year")])
pred_glm1_1se$resid <- pred_glm1_1se[,response] - pred_glm1_1se[,paste0(response, "_pred")]
pred_glm1_1se$extremeResid <- NA
pred_glm1_1se[pred_glm1_1se$resid > 70 | pred_glm1_1se$resid < -70,"extremeResid"] <- 1Maps of Observations, Predictions, and Residuals=
Observations across the temporal range of the dataset
pred_glm1 <- pred_glm1 %>%
mutate(resid = .[[response]] - .[[paste0(response,"_pred")]])
# rasterize
# get reference raster
test_rast <- rast("../../../Data_raw/dayMet/rawMonthlyData/orders/70e0da02b9d2d6e8faa8c97d211f3546/Daymet_Monthly_V4R1/data/daymet_v4_prcp_monttl_na_1980.tif") %>%
terra::aggregate(fact = 8, fun = "mean")## |---------|---------|---------|---------|========================================= ## add ecoregion boundaries (for our ecoregion level model)
regions <- sf::st_read(dsn = "../../../Data_raw/Level2Ecoregions/", layer = "NA_CEC_Eco_Level2") ## Reading layer `NA_CEC_Eco_Level2' from data source
## `/Users/astears/Documents/Dropbox_static/Work/NAU_USGS_postdoc/PED_vegClimModels/Data_raw/Level2Ecoregions' using driver `ESRI Shapefile'
## Simple feature collection with 2261 features and 8 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -4334052 ymin: -3313739 xmax: 3324076 ymax: 4267265
## Projected CRS: Sphere_ARC_INFO_Lambert_Azimuthal_Equal_Arearegions <- regions %>%
st_transform(crs = st_crs(test_rast)) %>%
st_make_valid() #%>%
#st_crop(st_bbox(test_rast))
#
# goodRegions_temp <- st_overlaps(y = cropped_states, x = regions, sparse = FALSE) %>%
# rowSums()
# goodRegions <- regions[goodRegions_temp != 0,]
ecoregionLU <- data.frame("NA_L1NAME" = sort(unique(regions$NA_L1NAME)),
"newRegion" = c(NA, "Forest", "dryShrubGrass",
"dryShrubGrass", "Forest", "dryShrubGrass",
"dryShrubGrass", "Forest", "Forest",
"dryShrubGrass", "Forest", "Forest",
"Forest", "Forest", "dryShrubGrass",
NA
))
goodRegions <- regions %>%
left_join(ecoregionLU)
mapRegions <- goodRegions %>%
filter(!is.na(newRegion)) %>%
group_by(newRegion) %>%
summarise(geometry = sf::st_union(geometry)) %>%
ungroup() %>%
st_simplify(dTolerance = 1000)
#mapview(mapRegions)
# rasterize data
plotObs <- pred_glm1 %>%
drop_na(paste(response)) %>%
#slice_sample(n = 5e4) %>%
terra::vect(geom = c("Long", "Lat")) %>%
terra::set.crs(crs(test_rast)) %>%
terra::rasterize(y = test_rast,
field = response,
fun = mean) #%>%
#terra::aggregate(fact = 2, fun = mean, na.rm = TRUE) %>%
#terra::crop(ext(-1950000, 1000000, -1800000, 1000000))
# get the extent of this particular raster, and crop it accordingly
tempExt <- crds(plotObs, na.rm = TRUE)
plotObs_2 <- plotObs %>%
crop(ext(min(tempExt[,1]), max(tempExt[,1]),
min(tempExt[,2]), max(tempExt[,2]))
)
# make figures
ggplot() +
geom_spatraster(data = plotObs_2) +
geom_sf(data = mapRegions, fill = NA, col = "rosybrown4", lwd = .5) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% st_crop(st_bbox(plotObs_2)),fill=NA ) +
labs(title = paste0("Observations of ", response, " in the ",ecoregion, " ecoregion")) +
scale_fill_gradient2(low = "brown",
mid = "wheat" ,
high = "darkgreen" ,
midpoint = 0, na.value = "lightgrey") +
xlim(st_bbox(plotObs_2)[c(1,3)]) +
ylim(st_bbox(plotObs_2)[c(2,4)])
Predictions across the temporal range of the dataset
# rasterize data
plotPred <- pred_glm1 %>%
drop_na(paste0(response,"_pred")) %>%
#slice_sample(n = 5e4) %>%
terra::vect(geom = c("Long", "Lat")) %>%
terra::set.crs(crs(test_rast)) %>%
terra::rasterize(y = test_rast,
field = paste0(response,"_pred"),
fun = mean) #%>%
#terra::aggregate(fact = 2, fun = mean, na.rm = TRUE) %>%
#terra::crop(ext(-1950000, 1000000, -1800000, 1000000))
# get the point location of those predictions that are > 100
highPred_points <- pred_glm1 %>%
filter(.[[paste0(response,"_pred")]] > 100 |
.[[paste0(response, "_pred")]] < 0) %>%
terra::vect(geom = c("Long", "Lat")) %>%
terra::set.crs(crs(test_rast))
# get the extent of this particular raster, and crop it accordingly
tempExt <- crds(plotPred, na.rm = TRUE)
plotPred_2 <- plotPred %>%
crop(ext(min(tempExt[,1]), max(tempExt[,1]),
min(tempExt[,2]), max(tempExt[,2]))
)
# make figures
ggplot() +
geom_spatraster(data = plotPred_2) +
geom_sf(data = mapRegions, fill = NA, col = "rosybrown4", lwd = .5) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% st_crop(st_bbox(plotObs_2)),fill=NA ) +
geom_sf(data = highPred_points, col = "red") +
labs(title = paste0("Predictions from the 'best lambda' fitted model of ", response, " in the ",ecoregion, " ecoregion"),
subtitle = "bestLambda model") +
scale_fill_gradient2(low = "wheat",
mid = "darkgreen",
high = "red" ,
midpoint = 100, na.value = "lightgrey",
limits = c(0,100)) +
xlim(st_bbox(plotObs_2)[c(1,3)]) +
ylim(st_bbox(plotObs_2)[c(2,4)])
# rasterize data
plotPred <- pred_glm1_1se %>%
drop_na(paste0(response,"_pred")) %>%
#slice_sample(n = 5e4) %>%
terra::vect(geom = c("Long", "Lat")) %>%
terra::set.crs(crs(test_rast)) %>%
terra::rasterize(y = test_rast,
field = paste0(response,"_pred"),
fun = mean) #%>%
#terra::aggregate(fact = 2, fun = mean, na.rm = TRUE) %>%
#terra::crop(ext(-1950000, 1000000, -1800000, 1000000))
# get the point location of those predictions that are > 100
highPred_points <- pred_glm1_1se %>%
filter(.[[paste0(response,"_pred")]] > 100 |
.[[paste0(response, "_pred")]] < 0) %>%
terra::vect(geom = c("Long", "Lat")) %>%
terra::set.crs(crs(test_rast))
# get the extent of this particular raster, and crop it accordingly
tempExt <- crds(plotPred, na.rm = TRUE)
plotPred_2 <- plotPred %>%
crop(ext(min(tempExt[,1]), max(tempExt[,1]),
min(tempExt[,2]), max(tempExt[,2]))
)
# make figures
ggplot() +
geom_spatraster(data = plotPred_2) +
geom_sf(data = mapRegions, fill = NA, col = "rosybrown4", lwd = .5) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% st_crop(st_bbox(plotObs_2)),fill=NA ) + geom_sf(data = highPred_points, col = "red") +
labs(title = paste0("Predictions from the '1SE lambda' fitted model of ", response, " in the ",ecoregion, " ecoregion"),
subtitle = "1 SE Lambda model") +
scale_fill_gradient2(low = "wheat",
mid = "darkgreen",
high = "red" ,
midpoint = 100, na.value = "lightgrey",
limits = c(0,100)) +
xlim(st_bbox(plotObs_2)[c(1,3)]) +
ylim(st_bbox(plotObs_2)[c(2,4)])
Residuals across the entire temporal extent of the dataset
# rasterize data
plotResid_rast <- pred_glm1 %>%
drop_na(resid) %>%
#slice_sample(n = 5e4) %>%
terra::vect(geom = c("Long", "Lat")) %>%
terra::set.crs(crs(test_rast)) %>%
terra::rasterize(y = test_rast,
field = "resid",
fun = mean) #%>%
#terra::aggregate(fact = 2, fun = mean, na.rm = TRUE) %>%
#terra::crop(ext(-1950000, 1000000, -1800000, 1000000))
# get the extent of this particular raster, and crop it accordingly
tempExt <- crds(plotResid_rast, na.rm = TRUE)
plotResid_rast_2 <- plotResid_rast %>%
crop(ext(min(tempExt[,1]), max(tempExt[,1]),
min(tempExt[,2]), max(tempExt[,2]))
)
# identify locations where residuals are >100 or < -100
badResids_high <- pred_glm1 %>%
filter(resid > 100) %>%
terra::vect(geom = c("Long", "Lat")) %>%
terra::set.crs(crs(test_rast))
badResids_low <- pred_glm1 %>%
filter(resid < -100) %>%
terra::vect(geom = c("Long", "Lat")) %>%
terra::set.crs(crs(test_rast))
# make figures
map <- ggplot() +
geom_spatraster(data =plotResid_rast_2) +
geom_sf(data = mapRegions, fill = NA, col = "rosybrown4", lwd = .5) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% st_crop(st_bbox(plotObs_2)),fill=NA ) +
geom_sf(data = badResids_high, col = "blue") +
geom_sf(data = badResids_low, col = "red") +
labs(title = paste0("Resids. (obs. - pred.) from Grass/shrub ecoregion-wide model of ", response),
subtitle = "bestLambda model \n red points indicate locations that have residuals below -100 \n blue points indicate locatiosn that have residuals above 100") +
scale_fill_gradient2(low = "red",
mid = "white" ,
high = "blue" ,
midpoint = 0, na.value = "lightgrey",
limits = c(-100,100)
) +
xlim(st_bbox(plotObs_2)[c(1,3)]) +
ylim(st_bbox(plotObs_2)[c(2,4)])
hist <- ggplot(pred_glm1) +
geom_histogram(aes(resid), fill = "lightgrey", col = "darkgrey") +
geom_text(aes(x = min(resid)*.9, y = 1500, label = paste0("min = ", round(min(resid),2)))) +
geom_text(aes(x = max(resid)*.9, y = 1500, label = paste0("max = ", round(max(resid),2))))
library(ggpubr)
ggarrange(map, hist, heights = c(3,1), ncol = 1)
# rasterize data
plotResid_rast <- pred_glm1_1se %>%
drop_na(resid) %>%
#slice_sample(n = 5e4) %>%
terra::vect(geom = c("Long", "Lat")) %>%
terra::set.crs(crs(test_rast)) %>%
terra::rasterize(y = test_rast,
field = "resid",
fun = mean) #%>%
#terra::aggregate(fact = 2, fun = mean, na.rm = TRUE) %>%
#terra::crop(ext(-1950000, 1000000, -1800000, 1000000))
# get the extent of this particular raster, and crop it accordingly
tempExt <- crds(plotResid_rast, na.rm = TRUE)
plotResid_rast_2 <- plotResid_rast %>%
crop(ext(min(tempExt[,1]), max(tempExt[,1]),
min(tempExt[,2]), max(tempExt[,2]))
)
# identify locations where residuals are >100 or < -100
badResids_high <- pred_glm1_1se %>%
filter(resid > 100) %>%
terra::vect(geom = c("Long", "Lat")) %>%
terra::set.crs(crs(test_rast))
badResids_low <- pred_glm1_1se %>%
filter(resid < -100) %>%
terra::vect(geom = c("Long", "Lat")) %>%
terra::set.crs(crs(test_rast))
# make figures
map <- ggplot() +
geom_spatraster(data =plotResid_rast_2) +
geom_sf(data = mapRegions, fill = NA, col = "rosybrown4", lwd = .5) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% st_crop(st_bbox(plotObs_2)),fill=NA ) +
geom_sf(data = badResids_high, col = "blue") +
geom_sf(data = badResids_low, col = "red") +
labs(title = paste0("Resids. (obs. - pred.) from Grass/shrub ecoregion-wide model of ", response),
subtitle = "1 SE Lambda model \n red points indicate locations that have residuals below -100 \n blue points indicate locatiosn that have residuals above 100") +
scale_fill_gradient2(low = "red",
mid = "white" ,
high = "blue" ,
midpoint = 0, na.value = "lightgrey",
limits = c(-100,100)
) +
xlim(st_bbox(plotObs_2)[c(1,3)]) +
ylim(st_bbox(plotObs_2)[c(2,4)])
hist <- ggplot(pred_glm1_1se) +
geom_histogram(aes(resid), fill = "lightgrey", col = "darkgrey") +
geom_text(aes(x = min(resid)*.9, y = 1500, label = paste0("min = ", round(min(resid),2)))) +
geom_text(aes(x = max(resid)*.9, y = 1500, label = paste0("max = ", round(max(resid),2))))
ggarrange(map, hist, heights = c(3,1), ncol = 1)
### Are there biases of the model predictions across year/lat/long?
# plot residuals against Year
yearResidMod_bestLambda <- ggplot(pred_glm1) +
geom_point(aes(x = jitter(Year), y = resid), alpha = .1) +
geom_smooth(aes(x = Year, y = resid)) +
xlab("Year") +
ylab("Residual from best lambda model") +
ggtitle("from best lamba model")
yearResidMod_1seLambda <- ggplot(pred_glm1_1se) +
geom_point(aes(x = jitter(Year), y = resid), alpha = .1) +
geom_smooth(aes(x = Year, y = resid)) +
xlab("Year") +
ylab("Residual from 1 SE lambda model")+
ggtitle("from 1 SE lamba model")
# plot residuals against Lat
latResidMod_bestLambda <- ggplot(pred_glm1) +
geom_point(aes(x = Lat, y = resid), alpha = .1) +
geom_smooth(aes(x = Lat, y = resid)) +
xlab("Latitude") +
ylab("Residual from best lambda model") +
ggtitle("from best lamba model")
latResidMod_1seLambda <- ggplot(pred_glm1_1se) +
geom_point(aes(x = Lat, y = resid), alpha = .1) +
geom_smooth(aes(x = Lat, y = resid)) +
xlab("Latitude") +
ylab("Residual from 1 SE lambda model") +
ggtitle("from 1 SE lamba model")
# plot residuals against Long
longResidMod_bestLambda <- ggplot(pred_glm1) +
geom_point(aes(x = Long, y = resid), alpha = .1) +
geom_smooth(aes(x = Long, y = resid)) +
xlab("Longitude") +
ylab("Residual from best lambda model") +
ggtitle("from best lamba model")
longResidMod_1seLambda <- ggplot(pred_glm1_1se) +
geom_point(aes(x = Long, y = resid), alpha = .1) +
geom_smooth(aes(x = Long, y = resid)) +
xlab("Longitude") +
ylab("Residual from 1 SE lambda model") +
ggtitle("from 1 SE lamba model")
library(patchwork)
(yearResidMod_bestLambda + yearResidMod_1seLambda) /
( latResidMod_bestLambda + latResidMod_1seLambda) /
( longResidMod_bestLambda + longResidMod_1seLambda)
Quantile plots
Binning predictor variables into “Deciles” (actually percentiles) and looking at the mean predicted probability for each percentile. The use of the word Decentiles is just a legacy thing (they started out being actual Percentiles)
# get deciles for best lambda model
if (length(prednames_fig) == 0) {
print("The best lambda model only contains one predictor (an intercept), so decile plots aren't possible to generate")
} else {
pred_glm1_deciles <- predvars2deciles(pred_glm1,
response_vars = response_vars,
pred_vars = prednames_fig,
cut_points = seq(0, 1, 0.005))
}
# get deciles for 1 SE lambda model
if (length(prednames_fig_1se) == 0) {
print("The 1SE lambda model only contains one predictor (an intercept), so decile plots aren't possible to generate")
} else {
pred_glm1_deciles_1se <- predvars2deciles(pred_glm1_1se,
response_vars = response_vars,
pred_vars = prednames_fig_1se,
cut_points = seq(0, 1, 0.005))
}Here is a quick version of LOESS curves fit to raw data (to double-check the quantile plot calculations)
if (length(prednames_fig) == 0) {
print("The model only contains one predictor (an intercept), so decile plots aren't possible to generate")
} else {
pred_glm1 %>%
select(all_of(c(prednames_fig, response_vars))) %>%
pivot_longer(cols = prednames_fig) %>%
ggplot() +
facet_wrap(~name, scales = "free") +
geom_point(aes(x = value, y = .data[[paste(response)]]), col = "darkblue", alpha = .1) + # observed values
geom_point(aes(x = value, y = .data[[response_vars[2]]]), col = "lightblue", alpha = .1) + # model-predicted values
geom_smooth(aes(x = value, y = .data[[paste(response)]]), col = "black", se = FALSE) +
geom_smooth(aes(x = value, y = .data[[response_vars[2]]]), col = "brown", se = FALSE)
}
Below are the actual quantile plots (note that the predictor variables are scaled)
if (length(prednames_fig) == 0) {
print("The model only contains one predictor (an intercept), so decile plots aren't possible to generate")
} else {
# publication quality version
g3 <- decile_dotplot_pq(pred_glm1_deciles, response = response, IQR = TRUE) + ggtitle("Decile Plot")
g4 <- add_dotplot_inset(g3, df = pred_glm1_deciles, dfRaw = pred_glm1, add_smooth = TRUE, deciles = FALSE)
if(save_figs) {
png(paste0("figures/quantile_plots/quantile_plot_", response, "_",ecoregion,".png"),
units = "in", res = 600, width = 5.5, height = 3.5 )
print(g4)
dev.off()
}
g4
}
if (length(prednames_fig_1se) == 0) {
print("The 1 se lambda model only contains one predictor (an intercept), so decile plots aren't possible to generate")
} else {
# publication quality version
g3 <- decile_dotplot_pq(pred_glm1_deciles_1se, response = response, IQR = TRUE) + ggtitle("Decile Plot")
g4 <- add_dotplot_inset(g3, df = pred_glm1_deciles_1se, dfRaw = pred_glm1_1se, add_smooth = TRUE, deciles = FALSE)
if(save_figs) {
png(paste0("figures/quantile_plots/quantile_plot_", response, "_",ecoregion,".png"),
units = "in", res = 600, width = 5.5, height = 3.5 )
print(g4)
dev.off()
}
g4
}
Deciles Filtered
20th and 80th percentiles for each climate variable
df <- pred_glm1[, prednames_fig] #%>%
#mutate(MAT = MAT - 273.15) # k to c
quantiles <- map(df, quantile, probs = c(0.2, 0.8), na.rm = TRUE)Filtered ‘Decile’ plots of data. These plots show each vegetation variable, but only based on data that falls into the upper and lower two deciles of each predictor variable.
if (length(prednames_fig) == 0) {
print("The model only contains one predictor (an intercept), so decile plots aren't possible to generate")
} else {
pred_glm1_deciles_filt <- predvars2deciles( pred_glm1,
response_vars = response_vars,
pred_vars = prednames_fig,
filter_var = TRUE,
filter_vars = prednames_fig,
cut_points = seq(0, 1, 0.005))
decile_dotplot_filtered_pq(pred_glm1_deciles_filt, xvars = prednames_fig)
#decile_dotplot_filtered_pq(pred_glm1_deciles_filt)
}## Processed 10976 groups out of 217178. 5% done. Time elapsed: 3s. ETA: 56s.Processed 16813 groups out of 217178. 8% done. Time elapsed: 4s. ETA: 47s.Processed 22649 groups out of 217178. 10% done. Time elapsed: 5s. ETA: 42s.Processed 28412 groups out of 217178. 13% done. Time elapsed: 6s. ETA: 39s.Processed 33864 groups out of 217178. 16% done. Time elapsed: 7s. ETA: 37s.Processed 39554 groups out of 217178. 18% done. Time elapsed: 8s. ETA: 35s.Processed 45379 groups out of 217178. 21% done. Time elapsed: 9s. ETA: 34s.Processed 51205 groups out of 217178. 24% done. Time elapsed: 10s. ETA: 32s.Processed 56804 groups out of 217178. 26% done. Time elapsed: 11s. ETA: 31s.Processed 62553 groups out of 217178. 29% done. Time elapsed: 12s. ETA: 29s.Processed 68157 groups out of 217178. 31% done. Time elapsed: 13s. ETA: 28s.Processed 73915 groups out of 217178. 34% done. Time elapsed: 14s. ETA: 27s.Processed 79754 groups out of 217178. 37% done. Time elapsed: 15s. ETA: 25s.Processed 85187 groups out of 217178. 39% done. Time elapsed: 16s. ETA: 24s.Processed 90803 groups out of 217178. 42% done. Time elapsed: 17s. ETA: 23s.Processed 96631 groups out of 217178. 44% done. Time elapsed: 18s. ETA: 22s.Processed 102476 groups out of 217178. 47% done. Time elapsed: 19s. ETA: 21s.Processed 108194 groups out of 217178. 50% done. Time elapsed: 20s. ETA: 20s.Processed 114027 groups out of 217178. 53% done. Time elapsed: 21s. ETA: 19s.Processed 119872 groups out of 217178. 55% done. Time elapsed: 22s. ETA: 17s.Processed 125708 groups out of 217178. 58% done. Time elapsed: 23s. ETA: 16s.Processed 131428 groups out of 217178. 61% done. Time elapsed: 24s. ETA: 15s.Processed 137221 groups out of 217178. 63% done. Time elapsed: 25s. ETA: 14s.Processed 142843 groups out of 217178. 66% done. Time elapsed: 26s. ETA: 13s.Processed 148592 groups out of 217178. 68% done. Time elapsed: 27s. ETA: 12s.Processed 154233 groups out of 217178. 71% done. Time elapsed: 28s. ETA: 11s.Processed 160065 groups out of 217178. 74% done. Time elapsed: 29s. ETA: 10s.Processed 165655 groups out of 217178. 76% done. Time elapsed: 30s. ETA: 9s.Processed 171090 groups out of 217178. 79% done. Time elapsed: 31s. ETA: 8s.Processed 176860 groups out of 217178. 81% done. Time elapsed: 32s. ETA: 7s.Processed 182436 groups out of 217178. 84% done. Time elapsed: 33s. ETA: 6s.Processed 188203 groups out of 217178. 87% done. Time elapsed: 34s. ETA: 5s.Processed 194050 groups out of 217178. 89% done. Time elapsed: 35s. ETA: 4s.Processed 199894 groups out of 217178. 92% done. Time elapsed: 36s. ETA: 3s.Processed 205646 groups out of 217178. 95% done. Time elapsed: 37s. ETA: 2s.Processed 211079 groups out of 217178. 97% done. Time elapsed: 38s. ETA: 1s.Processed 216921 groups out of 217178. 100% done. Time elapsed: 39s. ETA: 0s.Processed 217178 groups out of 217178. 100% done. Time elapsed: 39s. ETA: 0s.
Filtered quantile figure with middle 2 deciles also shown
if (length(prednames_fig) == 0) {
print("The model only contains one predictor (an intercept), so decile plots aren't possible to generate")
} else {
pred_glm1_deciles_filt_mid <- predvars2deciles(pred_glm1,
response_vars = response_vars,
pred_vars = prednames_fig,
filter_vars = prednames_fig,
filter_var = TRUE,
add_mid = TRUE,
cut_points = seq(0, 1, 0.005))
g <- decile_dotplot_filtered_pq(df = pred_glm1_deciles_filt_mid, xvars = prednames_fig)
g
if(save_figs) {x
jpeg(paste0("figures/quantile_plots/quantile_plot_filtered_mid_v1", , ".jpeg"),
units = "in", res = 600, width = 5.5, height = 6 )
g
dev.off()
}
}## Processed 15603 groups out of 325869. 5% done. Time elapsed: 3s. ETA: 59s.Processed 20942 groups out of 325869. 6% done. Time elapsed: 4s. ETA: 58s.Processed 25866 groups out of 325869. 8% done. Time elapsed: 5s. ETA: 57s.Processed 31188 groups out of 325869. 10% done. Time elapsed: 6s. ETA: 56s.Processed 36487 groups out of 325869. 11% done. Time elapsed: 7s. ETA: 55s.Processed 41896 groups out of 325869. 13% done. Time elapsed: 8s. ETA: 54s.Processed 47099 groups out of 325869. 14% done. Time elapsed: 9s. ETA: 53s.Processed 52157 groups out of 325869. 16% done. Time elapsed: 10s. ETA: 52s.Processed 57460 groups out of 325869. 18% done. Time elapsed: 11s. ETA: 51s.Processed 62840 groups out of 325869. 19% done. Time elapsed: 12s. ETA: 50s.Processed 68239 groups out of 325869. 21% done. Time elapsed: 13s. ETA: 49s.Processed 73445 groups out of 325869. 23% done. Time elapsed: 14s. ETA: 48s.Processed 78689 groups out of 325869. 24% done. Time elapsed: 15s. ETA: 47s.Processed 84080 groups out of 325869. 26% done. Time elapsed: 16s. ETA: 46s.Processed 89436 groups out of 325869. 27% done. Time elapsed: 17s. ETA: 45s.Processed 94873 groups out of 325869. 29% done. Time elapsed: 18s. ETA: 44s.Processed 100330 groups out of 325869. 31% done. Time elapsed: 19s. ETA: 42s.Processed 105776 groups out of 325869. 32% done. Time elapsed: 20s. ETA: 41s.Processed 111092 groups out of 325869. 34% done. Time elapsed: 21s. ETA: 40s.Processed 116484 groups out of 325869. 36% done. Time elapsed: 22s. ETA: 39s.Processed 121815 groups out of 325869. 37% done. Time elapsed: 23s. ETA: 38s.Processed 126810 groups out of 325869. 39% done. Time elapsed: 24s. ETA: 37s.Processed 131829 groups out of 325869. 40% done. Time elapsed: 25s. ETA: 37s.Processed 137229 groups out of 325869. 42% done. Time elapsed: 26s. ETA: 35s.Processed 142658 groups out of 325869. 44% done. Time elapsed: 27s. ETA: 34s.Processed 148100 groups out of 325869. 45% done. Time elapsed: 28s. ETA: 33s.Processed 153537 groups out of 325869. 47% done. Time elapsed: 29s. ETA: 32s.Processed 158907 groups out of 325869. 49% done. Time elapsed: 30s. ETA: 31s.Processed 164331 groups out of 325869. 50% done. Time elapsed: 31s. ETA: 30s.Processed 169660 groups out of 325869. 52% done. Time elapsed: 32s. ETA: 29s.Processed 174768 groups out of 325869. 54% done. Time elapsed: 33s. ETA: 28s.Processed 180203 groups out of 325869. 55% done. Time elapsed: 34s. ETA: 27s.Processed 184908 groups out of 325869. 57% done. Time elapsed: 35s. ETA: 27s.Processed 190519 groups out of 325869. 58% done. Time elapsed: 36s. ETA: 26s.Processed 196186 groups out of 325869. 60% done. Time elapsed: 37s. ETA: 25s.Processed 201770 groups out of 325869. 62% done. Time elapsed: 38s. ETA: 23s.Processed 207427 groups out of 325869. 64% done. Time elapsed: 39s. ETA: 22s.Processed 212805 groups out of 325869. 65% done. Time elapsed: 40s. ETA: 21s.Processed 218497 groups out of 325869. 67% done. Time elapsed: 41s. ETA: 20s.Processed 224212 groups out of 325869. 69% done. Time elapsed: 42s. ETA: 19s.Processed 229920 groups out of 325869. 71% done. Time elapsed: 43s. ETA: 18s.Processed 235597 groups out of 325869. 72% done. Time elapsed: 44s. ETA: 17s.Processed 240706 groups out of 325869. 74% done. Time elapsed: 45s. ETA: 16s.Processed 246193 groups out of 325869. 76% done. Time elapsed: 46s. ETA: 15s.Processed 251444 groups out of 325869. 77% done. Time elapsed: 47s. ETA: 14s.Processed 256824 groups out of 325869. 79% done. Time elapsed: 48s. ETA: 13s.Processed 262420 groups out of 325869. 81% done. Time elapsed: 49s. ETA: 12s.Processed 268047 groups out of 325869. 82% done. Time elapsed: 50s. ETA: 10s.Processed 273601 groups out of 325869. 84% done. Time elapsed: 51s. ETA: 9s.Processed 278782 groups out of 325869. 86% done. Time elapsed: 52s. ETA: 8s.Processed 284455 groups out of 325869. 87% done. Time elapsed: 53s. ETA: 7s.Processed 290130 groups out of 325869. 89% done. Time elapsed: 54s. ETA: 6s.Processed 295500 groups out of 325869. 91% done. Time elapsed: 55s. ETA: 5s.Processed 301184 groups out of 325869. 92% done. Time elapsed: 56s. ETA: 4s.Processed 306608 groups out of 325869. 94% done. Time elapsed: 57s. ETA: 3s.Processed 312287 groups out of 325869. 96% done. Time elapsed: 58s. ETA: 2s.Processed 317956 groups out of 325869. 98% done. Time elapsed: 59s. ETA: 1s.Processed 323509 groups out of 325869. 99% done. Time elapsed: 60s. ETA: 0s.Processed 325869 groups out of 325869. 100% done. Time elapsed: 61s. ETA: 0s.Cross-validation
Using best lambda model
Use terms from global model to re-fit and predict on different held out regions
Figures show residuals for each of the models fit to held-out ecoregions
These models were fit to six ecoregions, and then predict on the indicated heldout ecoregion
if (length(prednames_fig) == 0) {
print("The model only contains one predictor (an intercept), so cross validation isn't practical")
} else {
## code from Tredennick et al. 2020
# try each separate level II ecoregion as a test set
# make a list to hold output data
outList <- vector(mode = "list", length = length(sort(unique(modDat_1$NA_L2NAME))))
# obs_pred <- data.frame(ecoregion = character(),obs = numeric(),
# pred_opt = numeric(), pred_null = numeric()#,
# #pred_nopenalty = numeric()
# )
## get the model specification from the global model
mat <- as.matrix(coef(fit_glm_bestLambda, s = "lambda.min"))
mat2 <- mat[mat[,1] != 0,]
f_cv <- as.formula(paste0(response, " ~ ", paste0(names(mat2)[2:length(names(mat2))], collapse = " + ")))
X_cv <- model.matrix(object = f_cv, data = modDat_1_s)
# get response variable
y_cv <- as.matrix(modDat_1_s[,response])
# now, loop through so with each iteration, a different ecoregion is held out
for(i_eco in sort(unique(modDat_1_s$NA_L2NAME))){
# split into training and test sets
test_eco <- i_eco
print(test_eco)
# identify the rowID of observations to be in the training and test datasets
train <- which(modDat_1_s$NA_L2NAME!=test_eco) # data for all ecoregions that aren't 'i_eco'
test <- which(modDat_1_s$NA_L2NAME==test_eco) # data for the ecoregion that is 'i_eco'
trainDat_all <- modDat_1_s %>%
slice(train) %>%
select(-newRegion)
testDat_all <- modDat_1_s %>%
slice(test) %>%
select(-newRegion)
# get the model matrices for input and response variables for cross validation model specification
X_train <- as.matrix(X_cv[train,])
X_test <- as.matrix(X_cv[test,])
y_train <- modDat_1_s[train,response]
y_test <- modDat_1_s[test,response]
# get the model matrices for input and response variables for original model specification
X_train_glob <- as.matrix(X[train,])
X_test_glob <- as.matrix(X[test,])
y_train_glob <- modDat_1_s[train,response]
y_test_glob <- modDat_1_s[test,response]
train_eco <- modDat_1_s$NA_L2NAME[train]
## just try a regular glm w/ the components from the global model
fit_i <- glm(data = trainDat_all, formula = f_cv,
,
family = stats::Gamma(link = "log")
)
# lasso model predictions with the optimal lambda
optimal_pred <- predict(fit_i, newdata= testDat_all, type = "response"
)
# null model and predictions
# the "null" model in this case is the global model
# predict on the test data for this iteration w/ the global model
null_pred <- predict.glm(fit_glm_bestLambda, newdata = testDat_all, type = "response")
# save data
tmp <- data.frame(ecoRegion_holdout = rep(test_eco,length(y_test)),obs=y_test,
pred_opt=optimal_pred, pred_null=null_pred#,
#pred_nopenalty=nopen_pred
) %>%
cbind(testDat_all)
# calculate RMSE, bias, etc. of
# RMSE of CV model
RMSE_optimal <- yardstick::rmse(data = data.frame(optimal_pred, y_test), truth = "y_test", estimate = "optimal_pred")[1,]$.estimate
# RMSE of global model
RMSE_null <- yardstick::rmse(data = data.frame(null_pred, y_test), truth = "y_test", estimate = "null_pred")[1,]$.estimate
# bias of CV model
bias_optimal <- mean(y_test - optimal_pred)
# bias of global model
bias_null <- mean( y_test - null_pred )
# put output into a list
tmpList <- list("testRegion" = i_eco,
"modelObject" = fit_i,
"modelPredictions" = tmp,
"performanceMetrics" = data.frame("RMSE_cvModel" = RMSE_optimal,
"RMSE_globalModel" = RMSE_null,
"bias_cvModel" = bias_optimal,
"bias_globalModel" = bias_null))
# save model outputs
outList[[which(sort(unique(modDat_1_s$NA_L2NAME)) == i_eco)]] <- tmpList
}
}## [1] "COLD DESERTS"
## [1] "MEDITERRANEAN PIEDMONT"
## [1] "SEMIARID PLAIN AND PRAIRIES"
## [1] "TEMPERATE PRAIRIES"
## [1] "WARM DESERTS"
## [1] "WEST-CENTRAL SEMIARID PRAIRIES"Below are the RMSE and bias values for predictions made for each holdout level II ecoregion, compared to predictions from the global model for that same ecoregion
# table of model performance
map(outList, .f = function(x) {
cbind(data.frame("holdout region" = x$testRegion), x$performanceMetrics)
}
) %>%
purrr::list_rbind() %>%
kable(col.names = c("Held-out ecoregion", "RMSE of CV model", "RMSE of global model",
"bias of CV model - mean(obs-pred.)", "bias of global model- mean(obs-pred.)"),
caption = "Performance of Cross Validation using 'best lambda' model specification") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed")) | Held-out ecoregion | RMSE of CV model | RMSE of global model | bias of CV model - mean(obs-pred.) | bias of global model- mean(obs-pred.) |
|---|---|---|---|---|
| COLD DESERTS | 16.547761 | 14.619749 | 3.3204705 | 0.2943774 |
| MEDITERRANEAN PIEDMONT | 18.230763 | 13.833386 | -5.6663422 | -2.3936666 |
| SEMIARID PLAIN AND PRAIRIES | 9.206096 | 8.487211 | 2.3465454 | 0.6525862 |
| TEMPERATE PRAIRIES | 15.606439 | 2.450505 | -0.6749311 | 0.1935602 |
| WARM DESERTS | 17.435641 | 16.766420 | 2.2944775 | 0.3409402 |
| WEST-CENTRAL SEMIARID PRAIRIES | 9.761410 | 8.868124 | -4.1142164 | -0.8369391 |
# visualize model predictions
for (i in 1:length(unique(modDat_1_s$NA_L2NAME))) {
holdoutRegion <- outList[[i]]$testRegion
predictionData <- outList[[i]]$modelPredictions
modTerms <- as.matrix(coef(outList[[i]]$modelObject)) %>%
as.data.frame() %>%
filter(V1!=0) %>%
rownames()
# calculate residuals
predictionData <- predictionData %>%
mutate(resid = .[["obs"]] - .[["pred_opt"]] ,
resid_globMod = .[["obs"]] - .[["pred_null"]])
# rasterize
# use 'test_rast' from earlier
# rasterize data
plotObs <- predictionData %>%
drop_na(paste(response)) %>%
#slice_sample(n = 5e4) %>%
terra::vect(geom = c("Long", "Lat")) %>%
terra::set.crs(crs(test_rast)) %>%
terra::rasterize(y = test_rast,
field = "resid",
fun = mean) #%>%
#terra::aggregate(fact = 2, fun = mean, na.rm = TRUE) %>%
#terra::crop(ext(-1950000, 1000000, -1800000, 1000000))
tempExt <- crds(plotObs, na.rm = TRUE)
plotObs_2 <- plotObs %>%
crop(ext(min(tempExt[,1]), max(tempExt[,1]),
min(tempExt[,2]), max(tempExt[,2]))
)
# identify locations where residuals are >100 or < -100
badResids_high <- predictionData %>%
filter(resid > 100) %>%
terra::vect(geom = c("Long", "Lat")) %>%
terra::set.crs(crs(test_rast))
badResids_low <- predictionData %>%
filter(resid < -100) %>%
terra::vect(geom = c("Long", "Lat")) %>%
terra::set.crs(crs(test_rast))
# make figures
# make histogram
hist_i <- ggplot(predictionData) +
geom_histogram(aes(resid_globMod), col = "darkgrey", fill = "lightgrey") +
xlab(c("Residuals (obs. - pred.)"))
# make map
map_i <- ggplot() +
geom_spatraster(data = plotObs_2) +
geom_sf(data = mapRegions, fill = NA, col = "rosybrown4", lwd = .5) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% st_crop(st_bbox(plotObs_2)),fill=NA ) +
geom_sf(data = badResids_high, col = "blue") +
geom_sf(data = badResids_low, col = "red") +
labs(title = paste0("Residuals (obs. - pred.) for predictions of \n", holdoutRegion, " \n from a model fit to other ecoregions"),
subtitle = paste0(response, " ~ ", paste0( modTerms, collapse = " + "))) +
scale_fill_gradient2(low = "red",
mid = "white" ,
high = "blue" ,
midpoint = 0, na.value = "lightgrey",
limits = c(-100, 100)) +
xlim(st_bbox(plotObs_2)[c(1,3)]) +
ylim(st_bbox(plotObs_2)[c(2,4)])
assign(paste0("residPlot_",holdoutRegion),
value = ggarrange(map_i, hist_i, heights = c(3,1), ncol = 1)
)
}
lapply(unique(modDat_1_s$NA_L2NAME), FUN = function(x) {
get(paste0("residPlot_", x))
})## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
Using 1se lambda model
Use terms from global model to re-fit and predict on different held out regions
Figures show residuals for each of the models fit to held-out ecoregions
These models were fit to six ecoregions, and then predict on the indicated heldout ecoregion
if (length(prednames_fig_1se) == 0) {
print("The model only contains one predictor (an intercept), so cross validation isn't practical")
} else {
## code from Tredennick et al. 2020
# try each separate level II ecoregion as a test set
# make a list to hold output data
outList <- vector(mode = "list", length = length(sort(unique(modDat_1$NA_L2NAME))))
# obs_pred <- data.frame(ecoregion = character(),obs = numeric(),
# pred_opt = numeric(), pred_null = numeric()#,
# #pred_nopenalty = numeric()
# )
## get the model specification from the global model
mat <- as.matrix(coef(fit_glm_se, s = "lambda.min"))
mat2 <- mat[mat[,1] != 0,]
f_cv <- as.formula(paste0(response, " ~ ", paste0(names(mat2)[2:length(names(mat2))], collapse = " + ")))
X_cv <- model.matrix(object = f_cv, data = modDat_1_s)
# get response variable
y_cv <- as.matrix(modDat_1_s[,response])
# now, loop through so with each iteration, a different ecoregion is held out
for(i_eco in sort(unique(modDat_1_s$NA_L2NAME))){
# split into training and test sets
test_eco <- i_eco
print(test_eco)
# identify the rowID of observations to be in the training and test datasets
train <- which(modDat_1_s$NA_L2NAME!=test_eco) # data for all ecoregions that aren't 'i_eco'
test <- which(modDat_1_s$NA_L2NAME==test_eco) # data for the ecoregion that is 'i_eco'
trainDat_all <- modDat_1_s %>%
slice(train) %>%
select(-newRegion)
testDat_all <- modDat_1_s %>%
slice(test) %>%
select(-newRegion)
# get the model matrices for input and response variables for cross validation model specification
X_train <- as.matrix(X_cv[train,])
X_test <- as.matrix(X_cv[test,])
y_train <- modDat_1_s[train,response]
y_test <- modDat_1_s[test,response]
# get the model matrices for input and response variables for original model specification
X_train_glob <- as.matrix(X[train,])
X_test_glob <- as.matrix(X[test,])
y_train_glob <- modDat_1_s[train,response]
y_test_glob <- modDat_1_s[test,response]
train_eco <- modDat_1_s$NA_L2NAME[train]
## just try a regular glm w/ the components from the global model
fit_i <- glm(data = trainDat_all, formula = f_cv,
,
family = stats::Gamma(link = "log")
)
coef(fit_i)
# lasso model predictions with the optimal lambda
optimal_pred <- predict(fit_i, newdata= testDat_all, type = "response"
)
# null model and predictions
# the "null" model in this case is the global model
# predict on the test data for this iteration w/ the global model
null_pred <- predict.glm(fit_glm_se, newdata = testDat_all, type = "response")
# save data
tmp <- data.frame(ecoRegion_holdout = rep(test_eco,length(y_test)),obs=y_test,
pred_opt=optimal_pred, pred_null=null_pred#,
#pred_nopenalty=nopen_pred
) %>%
cbind(testDat_all)
# calculate RMSE, bias, etc. of
# RMSE of CV model
RMSE_optimal <- yardstick::rmse(data = data.frame(optimal_pred, y_test), truth = "y_test", estimate = "optimal_pred")[1,]$.estimate
# RMSE of global model
RMSE_null <- yardstick::rmse(data = data.frame(null_pred, y_test), truth = "y_test", estimate = "null_pred")[1,]$.estimate
# bias of CV model
bias_optimal <- mean(y_test - optimal_pred)
# bias of global model
bias_null <- mean(y_test - null_pred )
# put output into a list
tmpList <- list("testRegion" = i_eco,
"modelObject" = fit_i,
"modelPredictions" = tmp,
"performanceMetrics" = data.frame("RMSE_cvModel" = RMSE_optimal,
"RMSE_globalModel" = RMSE_null,
"bias_cvModel" = bias_optimal,
"bias_globalModel" = bias_null))
# save model outputs
outList[[which(sort(unique(modDat_1_s$NA_L2NAME)) == i_eco)]] <- tmpList
}
}## [1] "COLD DESERTS"
## [1] "MEDITERRANEAN PIEDMONT"
## [1] "SEMIARID PLAIN AND PRAIRIES"
## [1] "TEMPERATE PRAIRIES"
## [1] "WARM DESERTS"
## [1] "WEST-CENTRAL SEMIARID PRAIRIES"Below are the RMSE and bias values for predictions made for each holdout level II ecoregion, compared to predictions from the global model for that same ecoregion
if (length(prednames_fig_1se) == 0) {
print("The model only contains one predictor (an intercept), so cross validation isn't practical")
} else {
# table of model performance
map(outList, .f = function(x) {
cbind(data.frame("holdout region" = x$testRegion), x$performanceMetrics)
}
) %>%
purrr::list_rbind() %>%
kable(col.names = c("Held-out ecoregion", "RMSE of CV model", "RMSE of global model",
"bias of CV model - mean(obs-pred.)", "bias of global model - mean(obs-pred.)"),
caption = "Performance of Cross Validation using '1 SE lambda' model specification") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
}| Held-out ecoregion | RMSE of CV model | RMSE of global model | bias of CV model - mean(obs-pred.) | bias of global model - mean(obs-pred.) |
|---|---|---|---|---|
| COLD DESERTS | 16.289655 | 15.179813 | 4.9546701 | 0.6624492 |
| MEDITERRANEAN PIEDMONT | 15.623357 | 15.188536 | -5.5078580 | -4.5737635 |
| SEMIARID PLAIN AND PRAIRIES | 9.151986 | 8.813295 | 0.8376364 | 0.6200745 |
| TEMPERATE PRAIRIES | 3.567955 | 2.401475 | 0.5737164 | 0.3600259 |
| WARM DESERTS | 17.955609 | 17.665875 | -2.5610153 | -1.3603285 |
| WEST-CENTRAL SEMIARID PRAIRIES | 9.301570 | 9.007804 | -2.4475253 | -0.7824258 |
if (length(prednames_fig_1se) == 0) {
print("The model only contains one predictor (an intercept), so cross validation isn't practical")
} else {
for (i in 1:length(unique(modDat_1_s$NA_L2NAME))) {
holdoutRegion <- outList[[i]]$testRegion
predictionData <- outList[[i]]$modelPredictions
modTerms <- as.matrix(coef(outList[[i]]$modelObject)) %>%
as.data.frame() %>%
filter(V1!=0) %>%
rownames()
# calculate residuals
predictionData <- predictionData %>%
mutate(resid = .[["obs"]] - .[["pred_opt"]] ,
resid_globMod = .[["obs"]] - .[["pred_null"]])
# rasterize
# use 'test_rast' from earlier
# rasterize data
plotObs <- predictionData %>%
drop_na(paste(response)) %>%
#slice_sample(n = 5e4) %>%
terra::vect(geom = c("Long", "Lat")) %>%
terra::set.crs(crs(test_rast)) %>%
terra::rasterize(y = test_rast,
field = "resid",
fun = mean) #%>%
#terra::aggregate(fact = 2, fun = mean, na.rm = TRUE) %>%
#terra::crop(ext(-1950000, 1000000, -1800000, 1000000))
tempExt <- crds(plotObs, na.rm = TRUE)
plotObs_2 <- plotObs %>%
crop(ext(min(tempExt[,1]), max(tempExt[,1]),
min(tempExt[,2]), max(tempExt[,2]))
)
# identify locations where residuals are >100 or < -100
badResids_high <- predictionData %>%
filter(resid > 100) %>%
terra::vect(geom = c("Long", "Lat")) %>%
terra::set.crs(crs(test_rast))
badResids_low <- predictionData %>%
filter(resid < -100) %>%
terra::vect(geom = c("Long", "Lat")) %>%
terra::set.crs(crs(test_rast))
# make figures
# make histogram
hist_i <- ggplot(predictionData) +
geom_histogram(aes(resid_globMod), col = "darkgrey", fill = "lightgrey") +
xlab(c("Residuals (obs. - pred.)"))
# make map
map_i <- ggplot() +
geom_spatraster(data = plotObs_2) +
geom_sf(data = mapRegions, fill = NA, col = "rosybrown4", lwd = .5) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% st_crop(st_bbox(plotObs_2)),fill=NA ) +
geom_sf(data = badResids_high, col = "blue") +
geom_sf(data = badResids_low, col = "red") +
labs(title = paste0("Residuals (obs. - pred.) for predictions of \n", holdoutRegion, " \n from a model fit to other ecoregions"),
subtitle = paste0(response, " ~ ", paste0( modTerms, collapse = " + "))) +
scale_fill_gradient2(low = "red",
mid = "white" ,
high = "blue" ,
midpoint = 0, na.value = "lightgrey",
limits = c(-100, 100)) +
xlim(st_bbox(plotObs_2)[c(1,3)]) +
ylim(st_bbox(plotObs_2)[c(2,4)])
assign(paste0("residPlot_",holdoutRegion),
value = ggarrange(map_i, hist_i, heights = c(3,1), ncol = 1)
)
}
lapply(unique(modDat_1_s$NA_L2NAME), FUN = function(x) {
get(paste0("residPlot_", x))
})
}## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
Save output
# # glm models
# #mods2save <- butcher::butcher(mod_glmFinal) # removes some model components so the saved object isn't huge
#
# #mods2save$formula <- best_form
# #mods2save$pred_vars_inter <- pred_vars_inter # so have interactions
# #n <- nrow(df_sample)
# #mods2save$data_rows <- n
#
#
# if(!test_run) {
# saveRDS(mods2save,
# paste0("./models/glm_beta_model_CONUSwide_", s, "_n", n,
# #sample_group,
# ".RDS"))
# if (byRegion == TRUE) {
# ## western forests
# saveRDS(mods2save_WF,
# paste0("./models/glm_beta_model_WesternForests_", s, "_n", n,
# #sample_group,
# ".RDS"))
# ## eastern forests
# saveRDS(mods2save_EF,
# paste0("./models/glm_beta_model_EasternForests_", s, "_n", n,
# #sample_group,
# ".RDS"))
# ## grass/shrub
# saveRDS(mods2save_G,
# paste0("./models/glm_beta_model_GrassShrub_", s, "_n", n,
# #sample_group,
# ".RDS"))
# }
# }## partial dependence plots
#vip::vip(mod_glmFinal, num_features = 15)
#pdp_all_vars(mod_glmFinal, mod_vars = pred_vars, ylab = 'probability',train = df_small)
#caret::varImp(fit)session info
Hash of current commit (i.e. to ID the version of the code used)
system("git rev-parse HEAD", intern=TRUE)## [1] "79890c55a196d40eb16ae968701c4515b44c260c"Packages etc.
sessionInfo()## R version 4.4.0 (2024-04-24)
## Platform: aarch64-apple-darwin20
## Running under: macOS Sonoma 14.7.5
##
## Matrix products: default
## BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: America/Denver
## tzcode source: internal
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] doMC_1.3.8 iterators_1.0.14 foreach_1.5.2 ggpubr_0.6.0 factoextra_1.0.7
## [6] USA.state.boundaries_1.0.1 glmnet_4.1-8 Matrix_1.7-0 kableExtra_1.4.0 rsample_1.2.1
## [11] here_1.0.1 StepBeta_2.1.0 ggtext_0.1.2 knitr_1.49 gridExtra_2.3
## [16] pdp_0.8.2 GGally_2.2.1 lubridate_1.9.4 forcats_1.0.0 stringr_1.5.1
## [21] dplyr_1.1.4 purrr_1.0.4 readr_2.1.5 tidyr_1.3.1 tibble_3.2.1
## [26] tidyverse_2.0.0 caret_6.0-94 lattice_0.22-6 ggplot2_3.5.1 sf_1.0-20
## [31] tidyterra_0.6.1 terra_1.8-21 ggspatial_1.1.9 dtplyr_1.3.1 patchwork_1.3.0
##
## loaded via a namespace (and not attached):
## [1] RColorBrewer_1.1-3 rstudioapi_0.17.1 jsonlite_1.9.1 shape_1.4.6.1 magrittr_2.0.3 modeltools_0.2-23
## [7] farver_2.1.2 rmarkdown_2.29 vctrs_0.6.5 rstatix_0.7.2 htmltools_0.5.8.1 broom_1.0.7
## [13] Formula_1.2-5 pROC_1.18.5 sass_0.4.9 parallelly_1.37.1 KernSmooth_2.23-22 bslib_0.9.0
## [19] plyr_1.8.9 sandwich_3.1-0 zoo_1.8-12 cachem_1.1.0 commonmark_1.9.1 lifecycle_1.0.4
## [25] pkgconfig_2.0.3 R6_2.6.1 fastmap_1.2.0 future_1.33.2 digest_0.6.37 colorspace_2.1-1
## [31] furrr_0.3.1 rprojroot_2.0.4 pkgload_1.3.4 labeling_0.4.3 yardstick_1.3.1 timechange_0.3.0
## [37] mgcv_1.9-1 abind_1.4-8 compiler_4.4.0 proxy_0.4-27 aod_1.3.3 withr_3.0.2
## [43] backports_1.5.0 carData_3.0-5 betareg_3.1-4 DBI_1.2.3 ggstats_0.9.0 ggsignif_0.6.4
## [49] MASS_7.3-60.2 lava_1.8.0 classInt_0.4-10 gtools_3.9.5 ModelMetrics_1.2.2.2 tools_4.4.0
## [55] units_0.8-5 lmtest_0.9-40 future.apply_1.11.2 nnet_7.3-19 glue_1.8.0 nlme_3.1-164
## [61] gridtext_0.1.5 grid_4.4.0 reshape2_1.4.4 generics_0.1.3 recipes_1.1.0 gtable_0.3.6
## [67] tzdb_0.4.0 class_7.3-22 data.table_1.17.0 hms_1.1.3 utf8_1.2.4 car_3.1-2
## [73] xml2_1.3.7 flexmix_2.3-19 markdown_1.13 ggrepel_0.9.5 pillar_1.10.1 splines_4.4.0
## [79] survival_3.5-8 tidyselect_1.2.1 svglite_2.1.3 stats4_4.4.0 xfun_0.51 hardhat_1.4.0
## [85] timeDate_4032.109 stringi_1.8.4 yaml_2.3.10 evaluate_1.0.3 codetools_0.2-20 cli_3.6.4
## [91] rpart_4.1.23 systemfonts_1.2.1 munsell_0.5.1 jquerylib_0.1.4 Rcpp_1.0.14 globals_0.16.3
## [97] gower_1.0.1 listenv_0.9.1 viridisLite_0.4.2 ipred_0.9-15 scales_1.3.0 prodlim_2024.06.25
## [103] e1071_1.7-14 crayon_1.5.3 combinat_0.0-8 rlang_1.1.5 cowplot_1.1.3